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Abstract 

We present a collective coordinate approach to describe coupled phase oscil¬ 
lators. We apply the method to study synchronisation in a Kuramoto model. In 
our approach an A-dimensional Kuramoto model is reduced to an n-dimensional 
ordinary differential equation with n N, constituting an immense reduction 
in complexity. The onset of both local and global synchronisation is repro¬ 
duced to good numerical accuracy, and we are able to describe both soft and 
hard transitions. By introducing 2 collective coordinates the approach is able 
to describe the interaction of two partially synchronised clusters in the case of 
bimodally distributed native frequencies. Furthermore, our approach allows us 
to accurately describe finite size scalings of the critical coupling strength. We 
corroborate our analytical results by comparing with numerical simulations of 
the Kuramoto model with all-to-all coupling networks for several distributions 
of the native frequencies. 

Pacs numbers: 05.45.Xt, 89.75.-k, 89.75.Fb 


Despite their inherent complexity large networks of interacting dynam¬ 
ical entities often exhibit coordinated ordered behaviour such as mutual 
synchronisation. The macroscopic behaviour of complex networks arises 
as a complicated interplay between the dynamics of each microscopic node 
and the overall topological properties of the network. It is a formidable 
challenge to reduce the dynamics of large networks to a small number 
of active degrees of freedom that is capable of capturing these complex 
dynamical phenomena. This work is a contribution towards this goal. 

*georg.gottwald@sydney.edu.au 
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1 Introduction 


The collective behaviour of interacting oscillators in complex networks is ubiquitous 
in nature and has occupied scientists from as disparate areas as bi olog y, engineering, 
mathematics, physics and sociology for many years now [fsTlij]. These sys¬ 

tems often exhibit collective synchronisation whereby some or all oscillatory agents 
assume the same phase. Synchronisation behaviour is strongly dependent, amongst 
other factors, on the nature of the distribution of the native frequencies. In the case 
where all oscillators are connected with each other and where their native frequencies 
are unimodally distributed, for example, the onset of synchronisation as a function 
of the coupling strength is a soft transition, where the order parameter increases 
smoothly from zero as in a second-order phase transition. On the other hand, in 
the case of uniformly distributed frequencies, the onset of synchronisation is a hard 
transition, where at the critical coupling strength the order parameter has a non-zero 
value as in hrst-order phase transitions, with possible hysteresis (lol. 20 . S[n|. Cap¬ 
turing all these different dynamic behaviours is a challenging task. 


The collective behaviour of coupled oscillators such as synchronisation behaviour 
suggests that the dynamics of complex systems may (at least in certain cases) be 
described by a low dimensional dynamical system. To hnd these dimension-reduced 
descriptions is a formidable challenge with some remarkable results in recent years 
B 0, B, 0, 0]- work we propose a new approach to describe coupled 

phase oscillators and their non-trivial dynamics. Our approach is not restricted to 
a thermodynamic limit of inhnite many oscillators and allows for the study of hnite 
size effects B, BB 0 |. apparent in any real world networks. 


The particular approach proposed in this work seeks to hnd an approximate 
parametrisation of the s ync hronisation manifold by means of appropriately chosen 
collective coordinates (3, 14, [l5[ 16, B- The underlying premise is that the actual 


solution of the dynamical system assumes a specihc functional form the parameters 
of which are coined collective coordinates. The temporal evolution of the actual solu¬ 
tion is then described by the temporal evolution of those parameters, constituting an 
immense reduction in dimensionality. The functional form of the actual solution and 
the associated collective coordinates have to be specihed upon inspection of numerical 
simulations of the underlying system. For the Kuramoto model we will establish that 
the phases are linearly correlated with the native frequencies and we dehne the collec¬ 
tive coordinate to be the parameter relating the two. The method deals directly with 
the dynamical system rather than its associated macroscopic (inhnite-dimensional) 


description for the distribution or moments thereof l22|, ll3|, ll2|, l23|. It is non- 


perturbative in the sense that the solution is not written as an expansion in some 
small parameter. The paper is organized as follows. In Section [2] we introduce the 
Kuramoto model which constitutes a paradigm for studying coupled phase oscillators. 
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Our approach to achieve effective model reduction of the dynamics is introduced in 
Section El In Section 0] the method is applied to the Kuramoto model with all-to-all 
coupling with three different distributions for the native frequencies and we compare 
the results of direct numerical simulations of the full Kuramoto model with those 
of the proposed l-(or 2-)dimensional reduced model. We consider here a uniform 
native frequency distribution where a hard onset of synchronisation is experienced, a 
unimodal normal frequency distribution where a soft onset of synchronisation is ex¬ 
perienced, and thirdly a bimodal frequency distribution where global synchronisation 
is preceded by partial synchronisation of weakly coupled synchronised communities. 
We conclude with a summary and discussion in Section [5l 


2 Kuramoto model 


Weakly coupled limit cycle oscillators can be described in terms of their phases as an 
autonomous dynamical system. A widely used model which governs the dynamics of 
the phases ^ of N oscillators with native frequencies Ui is the celebrated Kuramoto 
model [13, 26, [ij 


K ^ 

^ ^ a-ij sin(9?j - ipi). 
i=i 


( 2 , 1 ) 


The adjacency matrix A = {a^} determines the topology of the network and describes 
which oscillators are connected. We restrict our analysis to unweighted, undirected 
networks for which the adjacency matrix A = {aij} is symmetric with Ojj = aji = 1 
if there is an edge between oscillators i and j, and a^j = 0 otherwise. The degree of a 
node di, i.e. the number of edges emanating from node i, is then given by di = a^j. 

For interacting oscillators, generically there exists a critical coupling strength Kc 
such that for sufficiently large coupling strength K > Kc the oscillators synchronise in 
the sense that they become locked to their mutual mean frequency and their phases 
become localized about their mean phase (l3. 18 . 26 |. This type of synchronous 


behaviour known as global synchronisation occurs if the dynamics settles on a globally 
attracting manifold p]. The level of synchronisation is often characterised by the order 
parameter [l3] 


N 


’'W = T7lE 




N' 


( 2 , 2 ) 


i=i 


with 0 < r < 1. In practice, the asymptotic limit of this order parameter 

^ r-To+T 

f = lim — / rit) dt , 

T-^-cx) T 


(2.3) 
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is estimated whereby Tq is chosen sufhciently large to eliminate transient behaviour 
of the oscillators. 

In the case of full synchronisation with (pi{t) = (pj{t) for all pairs i,j and for all 
times t we obtain f = r = 1. In the case where all oscillators behave independently 
with random initial conditions f = 0(1/\/N) indicates incoherent phase dynamics; 
values inbetween indicate partial coherence. 


3 Collective coordinate approach 


We will employ a non-perturbative approach to study synchronisation. Our approach 
is borrowed from the theory of solitary waves where it is known as collective coor¬ 
dinate approach 2^]: it has since been used in the context of dissipative pattern 
forming systems [3, 14, 15, [igI, 1^. The method we propose makes explicit use of the 
functional form of the phases as suggested by numerical simulations. The parameters 
describing the functional form of the phases constitute the collective coordinates. For 
example, if observations reveal that the functional form of the solution is bell-shaped 
at all times, the collective coordinates might be the amplitude and width of a Gaus¬ 
sian. The temporal evolution of the full solution is then described by the temporal 
evolution of the collective coordinates, i.e. how the amplitude and the width of the 
Gaussians evolve in time. Of course, a specihc assumed functional form is typically 
only an approximation of the actual solution. To eek out most of the assumed ansatz 
the collective coordinates are determined to optimally describe the solution. The 
most appropriate notion of optimality is to require that the error made by restricting 
the solution to be of the assumedz ansatz is minimised. Minimisation is achieved if 
the error is orthogonal to the subspace of the solutions spanned by the collective co¬ 
ordinates. This projection yields an evolution equation for the collective coordinates 
which allows to describe the actual solution at all times. 


We now establish the method of collective coordinates for the Kuramoto model 
in detail. Without loss of generality we assume that the mean frequency is zero 
(unless stated otherwise). Let us assume that the nodes are labelled in order of 
increasing native frequencies, i.e. i = 1 denotes the node with the most negative native 
frequency cui and i = N denotes the node with the most positive native frequency 
ujn- In Figure [U we show a snapshot of the phases (pj obtained by a numerical 
simulation of the Kuramoto model with an underlying Erdos-Renyi topology with 
N = 200 oscillators at a coupling strength K = 9.5. The associated order parameter 
is f = 0.78 indicating a high level of synchronisation. The hgure shows that the phases 
of oscillators with native frequencies of sufficiently small absolute value are frequency 
locked and correlate highly with the underlying native frequency distribution. This 
observation suggests that the phases of those frequency-locked oscillators may be 
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approximated by 


ipi{t) = a{t) cji . (3.1) 

Oscillators with large absolute native frequencies which could not be entrained at a 
given coupling strength, do not obey this functional relationship but rather oscillate 
with their native frequencies. The ansatz fl3.ip is trivially exact for K = 0 with 
a{t) = t. Furthermore, in the case of an all-to-all coupling the ansatz fl3.ip can be 
formally motivated for large coupling strength as follows. The stationary Kuramoto 
model fl2.ip can be rewritten as cuj = —Kr — ipi) with -0 being the mean phase 
[l^ . Expanding ipi = ^|J + arcsin(a;j/(riF)) in 1/K for large coupling strength yields 
up to hrst order ^pi = 'ip + uJi/{rK). Since the Kuramoto model fl2.1l) is invariant 
under constant phase shifts we may set 0 = 0 leading to our ansatz fl3.1l) . 

Our method consists of assuming that the phases of the N oscillators are approxi¬ 
mately given by our ansatz fl3.ip . The time-dependent amplitude a(t) takes the role 
of a collective coordinate. Our goal is to find an evolution equation for a{t) and 
thereby reducing the iV-dimensional Kuramoto model of phase oscillators to a one¬ 
dimensional ordinary differential equation for a{t) (in Section IT73] we will see how to 
modify the approach to include more collective coordinates). We do so by requiring 
that the error 


Sn — cnn. 


Ui 


K ^ 

-Y 

N ^ 


ttij sin(a(a;j — ooi)) 




made by restricting the solution to the subspace dehned by the ansatz (13.ip is min¬ 
imised. This is achieved by assuring that the error Sa is orthogonal to the restricted 
subspace spanned by fl3.ip . We therefore require that the error Sa is orthogonal to 
the tangent space of the solution manifold fl3.ip which is spanned by dipi/da = cuj. 
Projecting the error onto the restricted subspace spanned by (13.ip yields the desired 
evolution equation for a 


ct — 1 -|- 


K 1 


N N 

sin(a(a;j - ca,)) , 

i=l j=l 


(3.2) 


with 


1 ^ 

(3-3) 

i=i 

Solutions a* solving fl3.2p with d = hi correspond to phase-locked solutions rotat¬ 
ing uniformly with frequency hi and phases pj = a*Uj -|- Qt. The existence of such 
solutions corresponds to a synchronised state. The advantage of this approach is 
that it allows to study the onset of synchronisation of the iV-dimensional network by 
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analysing a one-dimensional problem and furthermore that it allows to study syn¬ 
chronisation for hnite network size N. 


In the limit N ^ oo we can simplify the expressions by introducing the frequency 
distribution g{uj) and the variance of the frequencies cx^ = hm 7 v-j>oo We obtain in 
an all-to-all coupling network with ajj = 1 for all i,j 


d — 1 -l- 


K 


(T. 


ug{u) 


sin(a !(?7 — uj))g{r})drjduj . 


(3.4) 


The order parameter r restricted to solutions ^Pj{t) = a{t)ujj is introduced as 


N 


= 


N ^ 


t=i 


(3.5) 


In the limit N ^ oo the real part yields 


f = 


cos{au)g{u) du, 


(3.6) 


where we used that our ansatz fl3.ip implies for the mean phase = 0. 

We remark that this approach is not restricted to all-to-all network topologies. 
For example, in an Erdos-Renyi network, where nodes are connected independently 
with probability p and where degrees dj are Poisson-distributed with mean degree 
d = pN, the inner sum in fl3.2p can be evaluated as a sum of (on average) d random 
variables pj ~ fi'(cu) with 

lim 

N^OO 




ttij sm 


'a{uj - Ui)) = d sm{a{p - uji))g{p)dp . 


The evolution equation for a{t) is then evaluated in the limit N ^ oo a.s 

ujg{uj) / sm{a{p — uj))g{p)dpdu . 


K 

d = 1 +p— 


(Tf. 


(3.7) 


In the next Section we will employ our framework to study the synchronisation 
properties of all-to-all coupling networks for several frequency distributions g{u). 


4 Examples 

We now set out to illustrate the capabilities of the collective coordinate approach to 
describe the synchronisation behaviour of phase oscillators in a Kuramoto model with 
an all-to-all coupling topology. We do so by determining the steady state solution a 
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Figure 1: Snapshot of the phases ipj for K = 9.5 for an = 200 obtained for the 
Kuramoto model with an Erdos-Renyi topology (with two nodes being connected 
with probability p = 0.05) and native frequencies oji = — 0.3 with ^ 

i 

The continuous line depicts a smooth cubic function. The corresponding value of the 
order parameter is f = 0.78. 

and the order parameter f, in the case of hnite N as well as in the thermodynamic 
limit of iV —)■ cx), for three different distributions of the native frequencies: uniform 
distribution, normal distribution and bimodal distribution. The results from the 
collective coordinate approach are then compared with results from direct numerical 
simulations of the corresponding Kuramoto model fl2.ip . 

Rather than performing averages over realisations of the native frequencies accord¬ 
ing to the respective distributions we will perform the calculations for the collective 
coordinate approach by choosing N values of the native frequencies such that the 
probability of a random draw of a native frequency to fall in the interval (cUjjCUj+i) is 
equal for all values of i. 

4.1 Uniform distribution of native frequencies 

In a hrst suite of experiments we consider native frequencies which are distributed 
uniformly on the interval [—1,1] with distribution 

g{uj) = 0.5 . (4.1) 

Dividing the compact support of the frequency distribution [—1,1] into iV —1 intervals 
of equal measure, i.e. Ui = 2{i — {N + l)/2)/(A^ — 1)) for i = 1, • • ■ , N, the evolution 
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equation fl3.2p for a for finite N is readily evaluated as 


a 


1 + -K - - - 

2 iV2(iV + l) 


CSC 


a 

N-1 


a ( 2Na 
N-1 iV- 1 


+ N sin 


2Na \ 
N-l) ’ 


(4.2) 


with 




N + 1 
3{N-1) ■ 


In the thermodynamic limit this simplihes to 

K sin a a cos a — sin a 

d = 1 H—r-- - - 

a 


(4.3) 


with al = limAT^oo = 1/3. 

The expression fl3.6p for the order parameter simplihes in the thermodynamic limit 
to 

sin(a) 

r =- 

a 

In Figure |2]we show the order parameter f as a function of the coupling strength K 
obtained from a long time integration of the full Kuramoto model (12.ip . The onset of 
synchronisation appears to be hard (see for example, Pazo 0) , i.e. there exists a non¬ 
zero value of the order parameter at the critical coupling strength Kc- The collective 
coordinate approach captures this very well as shown in Figure [2l Figure [3] shows 
that within the framework of collective coordinates the hard onset of synchronisation 
is described as a saddle node bifurcation j^: for K > Kc = 1.234 a pair of stationary 
solutions (fj = aujj (a smaller stable and a larger unstable one) exist; at criticality the 
two solutions collide in a saddle node bifurcation at a = ckc ~ 1.303, and there are 
no stationary solutions for K < Kc- Evaluating the right-hand-side of fl4.3p around 
the critical value ckc yields as an approximation of the stable and unstable stationary 
solutions CKs^u close to criticality 


(4.4) 


= Oc ± my 1 - ^ , (4-5) 

with the critical coupling strength Kc and m = 270/0.177. Figure [3] shows a 

numerical evaluation of the stationary solutions a of fl3.ip as well as the approximate 
solutions fl4.5p . Note that the stable stationary solution is well approximated for a 
large range of coupling strengths K even far away from criticality. 

We now analyse the order parameter r as given by fl4.4p . Figure [2] shows the or¬ 
der parameter as a function of the coupling strength obtained from a numerical 






















simulation of a large network with N = 10, 000 nodes simulating the Kuramoto 
model fl2.ip . and as calculated within the collective coordinate framework using fl4.4p . 
The critical coupling strengths for the full Kuramoto model with N = 10, 000 is 
Kc = 1.279 which is close to the exact analytical result for the thermodynamic limit 
with Tfc = d/vr ~ 1.273 (lol. 29 |. This is well approximated by our simple model 
with an error of 3%. The non-zero order parameter at the hard transition, which is 
Tc = Ti ~ 0.785 in the thermodynamic limit (lol. [^ . is estimated as = 0.744 
within the collective coordinate approach implying a 5% error. Note that the order 
parameter is extremely well approximated for large values of the coupling strength. 
This is not surprising since, as pointed out in Section [3l the collective coordinate 
ansatz fl3.ip is consistent with an expansion of the stationary solution in 1/77 for 
all-to-all coupling networks. 


A particular advantage of our approach is that it allows us to study the finite size 
scaling of synchronisation behaviour j^, l8[0, 27|. In Figure H] we show a comparison 
of the critical coupling strength Kc{N) as calculated via our collective coordinate 
approach for variable network sizes N and results from direct simulations of the Ku¬ 
ramoto model fl 2 .ip . The difficulty is determining the critical coupling strength Kc in 
hnite size networks is that the order parameter has fluctuations of the order 1 /\/N 
which confounds the onset. As a proxy for the critical coupling strength we record 
for each value of N the smallest value of the coupling strength K such that f > 0.8. 
We have also used the criterion whereby the critical coupling is determined as the 
coupling strength when the minimal value of the order parameter r{t) over some suf¬ 
ficiently long time window changes from values close to zero to values significantly 
above zero [28|. This method yields very pronounced onsets, but is not able to detect 
global synchronisation in the case when it is preceded by partial synchronisation. We 
therefore present only results obtained using the first method. In the case of uni¬ 
formly distributed native frequencies, however, both methods yield the same results. 
Linear regression suggests a scaling Kc{N) — K* ~ N where we estimate the critical 
coupling strength in the thermodynamic limit K* as K* = 1.279 for the full Ku¬ 
ramoto model and K* — 1.234 for the collective coordinate approach. 


Besides being able to describe the collective behaviour of oscillators and the onset 
of synchronisation, we now show that the collective coordinate approach also cap¬ 
tures the temporal evolution of individual oscillators through the evolution equation 
fl3.2p or its equivalent formulation fl4.2p for uniformly distributed native frequencies. 
For sufficiently small coupling strengths 77, where the oscillators only weakly inter¬ 
act, both models produce indistinguishable trajectories with phases growing linear in 
time (not shown). Figure E] shows a comparison of actual trajectories for a network 
with N = 101 oscillators at coupling strength 77 = 1.5 > Kc where the collective 
coordinate approach describes the order parameter f 0.9 very well (cf. Figure |2]). 
We show a comparison of the phase of the 75th oscillator 9975 with native frequency 
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Uj 5 = 0.48 is obtained by solving the full Kuranioto model fl2.ip and by solving fl4.2p 
for the collective coordinate approach fld.ip . If the initial conditions are chosen to 
satisfy (pj{0) = aoUj with the initial condition a(0) = ao not too far away from its 
equilibrium solution, the two trajectories are reasonably close (top panel). This cor¬ 
respondence of the time evolution of the solutions of the full Kuramoto model and 
the collective coordinate approach is destroyed for initial conditions which are too 
far from the asymptotic state, i.e. if ao is chosen too large. Their asymptotic state, 
however, will be close and both systems will evolve to the same hx point, implying 
that the order parameter f will be close for the two systems. Similarly, if the initial 
conditions (pj{0) of the Kuramoto model are distributed around the initial condition 
the asymptotic state and therefore differs from the initial condition (pj{0) implied 
by the collective coordinate ansatz fld.ip . the asymptotic temporal evolution of the 
full Kuramoto model and the reduced collective coordinate system are close (not 
shown). This is consistent with the previous observation that the order parameters 
f are close for the respective systems, as shown in Figure [2j We show a snapshot 
depicting the phases of all oscillators in the phase-locked state illustrating that the 
collective coordinate approach captures the dynamics of the full model. Deviations 
occur for the extreme oscillators with largest absolute value of the native frequencies. 
As we have seen in Figure |2] the collective coordinate approach predicts the onset 
of synchronisation for smaller values of K than observed for the actual Kuramoto 
model. For coupling strength where the order parameter signihcantly differs between 
the reduced model and the full model, there is of course, also no correspondence 
between the temporal evolution of the phases nor their asymptotic dynamics. We 
remark that we obtain similar results for networks differing in several orders of mag¬ 
nitude in size. For small networks of, for example, size N = 20, the phases are very 
well recovered if the native frequencies are chosen such that they divide the interval 
[—1,1] into equiprobable partitions. For a particular random draw from the uniform 
distribution, the phases and their asymptotic states may differ though, in particular 
for oscillators with large absolute native frequencies. This discrepancy can be allevi¬ 
ated for the well-synchronised oscillators if averages over many realisations of native 
frequencies are taken. With increasing size of the networks, the difference between 
solutions obtained for random realisations of the native frequencies become smaller. 


4.2 Normal distribution of native frequencies 

In a second suite of experiments, we consider native frequencies which are normally 
distributed with cjj ~ A/'(0, a^). The distribution is given by 

g(a.) = iexp(-^) , (4.6) 
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Figure 2: Order parameter f as a function of the coupling strength K for a network 
with uniformly distributed native frequencies. Depicted are results from a direct 
numerical integration of the Kuramoto model (12.11) with N = 10, 000 (crosses, online 
red) and from the collective coordinate approach fl4.4p (continuous line, online blue). 


with a normalisation constant Z = We use here = 0.1. The evolution 

equation fl3.2p for a for hnite N can be evaluated for random draws of cuj, but we omit 
here the cumbersome expressions. In the thermodynamic limit the dynamic model 
for the collective coordinate fl3.4p simplihes to 

a = 1 — Kaexp{—a^a^) . (4.7) 


The equation for the order parameter 03.61) can be evaluated in the thermodynamic 
limit to 


r = exp 


( 7 ,^, 0 !^ 


(4.8) 


It is well known that in the case of unimodal frequency distributions, the onset 


of synchronisation is soft jlOl. Il8| . This is illustrated in Figure [ 6 ] where f is shown as 
a function of the coupling strength. At the so called “Kuramoto coupling” K = Ki 
the order parameter becomes non-zero and a few oscillators with native frequencies 
close to the mean frequency mutually synchronise; increasing the coupling strength 
allows more and more oscillators to synchronise, implying a continuous change of the 
order parameter f{K) as supposed to the hard transition in the case of uniformly 
distributed native frequencies described in the previous subsection. At some coupling 
strength K = K^. global synchronisation sets in affecting all oscillators (3ol |. 

In the thermodynamic limit N ^ oo the Kuramoto coupling can be approximated 
by Ki = 2/Tig(0) 0.505 [l^. The transition to global synchronisation is not visible, 

however, by just looking at the order parameter f determined from numerical simu¬ 
lations of the full Kuramoto model fl2.ip . 

We will now show that the collective coordinate approach is able to describe both, 
the onset of global synchronisation at K = K,, as well as the onset of local synchro¬ 
nisation at the “Kuramoto coupling” K = Ki. The onset of global synchronisation 
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Figure 3: Top; Plot of the term ojf (dotted horizontal line, online red) and the 
term proportional to the coupling strength K in fl3.2l) as a function of the collective 
coordinate a for a network with uniformly distributed native frequencies. Intersec¬ 
tions denote stationary solutions of fl3.2p . Depicted are the subcritical case at iP = 1.1 
(dashed curve, online cyan), the critical case K = Kc = 1.234 (continuous curve, on¬ 
line blue) and the supercritical case with K = 1.4 (continuous line with circles, online 
magenta) for N = 1000. At criticality we hnd ac = 1.303. Bottom: The stable and 
unstable stationary solutions a as a function of the coupling strength K as calculated 
from the collective coordinate approach fl3.2p (continuous lines) and from the approx¬ 
imation fl4.5p (dashed lines). The two approximations are hardly distinguishable on 
the lower stable branch. 
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Figure 4: Scaling of the critical coupling strength Kc as a function of the network size 
N for a network with uniformly distributed native frequencies. Depicted are results 
from a direct numerical integration of the Kuramoto model fl2.ip (crosses, online red) 
and from the collective coordinate approach (continuous line, online blue). The two 
lines have slopes of 1. 




Figure 5: Phases <p(f) calculated from simulations of the full Kuramoto model fl2.ip 
(continuous lines, online red) and from the corresponding 1-dimensional system fl3.2l) 
for the collective coordinate with ipi = a{t)uji (crosses, online blue) for a network of 
N = 100 oscillators with uniformly distributed native frequencies at coupling strength 
K = 1.5. Top: Temporal evolution of for initial conditions (^,(0) = a^Ui with 

ao = 0.5. Bottom: Snapshot of the phases (pi{T) at time T = 20. 
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can be calculated as before. In Figure E] we show a result of the collective coordinate 
approach fl4.8p which predicts the onset of global synchronisation at k, 0.730 with 
a non-zero value of fc ~ 0.779. 

By construction, the ansatz fl3.ip cannot describe local synchronisation where only a 
subset of the N phase oscillators are phase locked. We now modify the collective co¬ 
ordinate approach to allow for local synchronisation. We denote by iV; the size of the 
mutually synchronised local group, consisting of those iV; oscillators with frequencies 
closest to the mean frequency zero. Hence we restrict our solutions to obey 

. ^ . s , N -Ni N + Ni , , 

= a{t)uj for -^- <j< -^-. (4.9) 

The evolution equation for the collective coordinate a{t) is again obtained by pro¬ 
jecting the error made by the ansatz fl4.9l) onto the restricted subspace spanned by 
04.91) . We obtain 


01 — 1 -|- 


K 1 


^12 ^12 

^ ^ sin(a(a;j - Ui)) , 

i=Nii j=Na 


where the variance of the local group of frequencies is 




1 


Ni2 

E 

j=i^ii 




(4.10) 


(4.11) 


with Nil = and W 2 = This is just the analogous formulation of 03.2p for 
a group of oscillators, centred around 0 ;^ = 0, of size iVj- Assuming that all those 
oscillators which can synchronise do so, the size of the locally synchronised group of 
oscillators Ni can be determined as the maximal value of Ni which supports stationary 
solutions of fl4.10p for a given coupling strength K. Note that Ni = N for K > Kc- 
Figure [7] shows how the normalised domain length of the local synchronised cluster 


L 


domain 


N ’ 


(4.12) 


increases from zero to ^domain > 0 ai K = Ki and then reaches ^domain = 1 at 
K = Kc at which point global synchronisation sets in. The Kuramoto coupling, i.e. 
the smallest value of K which gives rise to a non-zero value of Tdomaim is estimated 
for N = 1000 by our approach as Ki ^ 0.5 corresponding very well with the nu¬ 
merically observed onset of local synchronisation. The asymptotic value is given by 
Ki ^ 2/TTg{0) 0.505 [l^. It is pertinent to mention that in the case of uniformly 

distributed native frequencies, no stationary solutions a exist for any Ni < N, consis¬ 
tent with the absence of local synchronisation and the existence of a hard transition, 
as seen in Figure [2J 


14 

















Figure 6: Order parameter f as a function of the coupling strength K for a network 
with normally distributed native frequencies. Depicted are results from a direct nu¬ 
merical integration of the Kuramoto model fl2.ip with N = 1000 (continuous line, 
online red) and from the collective coordinate approach for global phase synchroni¬ 
sation (crosses, online blue) and for local phase synchronisation (open circles, online 
cyan). The curves coincide for sufficiently large values of the coupling strength K. 


In Figure [H] we illustrate again that the collective coordinate approach can be used 
to study hnite size scaling. For normally distributed native frequencies the numerics 
suggest a hnite size scaling of Kc{N) — K*{N) iV2/3. 

We show again a comparison of the actual temporal evolution of individual oscil¬ 
lators. Figure 0 shows results for the global synchronisation regime aX K = 0.9 and 
Figure dU] for the local synchronisation regime aX K = 0.6. In the case of the local 
synchronisation regime we assume that the oscillators which do not take part in the 
synchronised cluster are simply oscillating with their native frequencies and satisfy 
<pj(t) = (ao -l- t)ui. The temporal evolution is well described by the collective coordi¬ 
nate approach in both cases. It is clearly seen that, whereas the collective coordinate 
approach is able to capture the dynamics well of the well-entrained oscillators, it has 
difficulties describing the dynamics of the entrained extreme oscillators with large ab¬ 
solute native frequencies as seen in the insets of Figures |9] and [TOl This discrepancy is 
due to the collective coordinate approach, as employed here, not taking into account 
the interaction with the drifting extreme oscillators. 


4.3 Bimodal distribution of native frequencies 

In a third suite of experiments, we consider native frequencies which are distributed 
according to a bimodal distribution with maxima at m = ±D and 
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Figure 7: Normalised length of phase synchronised domain as a function of the cou¬ 
pling strength for a network with normal native frequency distribution, calculated 
using the collective coordinate approach. 



Figure 8: Scaling of the critical coupling strength as a function of the network size 
N for a network with normally distributed native frequencies. Depicted are results 
from a direct numerical integration of the Kuramoto model fl2.ip (crosses, online red) 
and from the collective coordinate approach (continuous line, online blue). The two 
lines have slopes of 0.664 suggesting a scaling with 2/3. 
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Figure 9: Phases ip{t) calculated from simulations of the full Kuramoto model 
fl2.1l) (continuous lines or open circles, online red) and from the corresponding 1- 
dimensional system fl3.2p for the collective coordinate with ipi = a{t)u)i (crosses, 
online blue) for a network of iV = 1000 oscillators with normally distributed native 
frequencies at coupling strength K = 0.9 corresponding to global synchronisation. 
Top: Temporal evolution of for initial conditions (^i(O) = aoUi with ao = 0.5. 

The native frequency is 07750 = 0.212. Bottom: Snapshot of the phases PiiT) at time 
T = 20. 




Figure 10: Phases p{t) calculated from simulations of the full Kuramoto model 
fl 2 .ll) (continuous lines or open circles, online red) and from the corresponding 1 - 
dimensional system fld.lOp for the collective coordinate with pi = a{t)ui (crosses or 
dashed line, online blue) for a network of = 1000 oscillators with normally dis¬ 
tributed native frequencies at coupling strength K = 0.6 with Ni = 805 correspond¬ 
ing to local synchronisation. Top: Temporal evolution of 97750 (t) for initial conditions 
(pj(0) = aoUi with ao = 0.5. The native frequency is 07750 = 0.212. Bottom: Snapshot 
of the phases Pi(T) at time T = 20. 
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Figure 11: Bimodal distribution g{u) of native frequencies fl4.13p with af, = 0.1 and 
n = 0.75. 


We choose here = 0.1 and = 0.75. The bimodal distribution for these parameters 
is depicted in Figure dH 

The synchronisation behaviour of Kuramoto networks with bimodal frequ ency 
distributions is more complex than in the two previous examples (lol. 0, Isl, 17, [l2|, 21 


If the two peaks are sufficiently close together, the behaviour is, roughly speaking, as 
described in the unimodal case, discussed in the previous section, with local synchro¬ 
nisation being organised by oscillators with native frequencies closest to the mean 
frequency zero. However, when the peaks are sufficiently separated, a so called stand¬ 
ing wave state occurs at some critical coupling strength K = Kp whereby the 
oscillators with native frequencies close to the peak frequencies ±11 may synchro¬ 
nise and form two synchronised clusters which rotate with the same frequency but in 
the opposite direction. Upon increasing the coupling strength further, the oscillators 
will eventually globally synchronise at a critical coupling strength K = Kc (lil. 21 


In Figure [12] we show a snapshot of the phases for the case Kp < K < Kc where 
two partially synchronised clusters are established, centred around the nodes with 
uji = ±11, respectively, which together form the standing wave state. In Figure [13] we 
show the order parameter f, where one can see clearly the standing wave state for 
Kp < K < Kc and global synchronisation for K > Kc with Kp ~ 1.05 and Kc ~ 1.7. 

First we apply our approach to the problem of global synchronisation, i.e. for 
K > 1.7. In the thermodynamic limit the dynamic model for the collective coordinate 
becomes 


(5 = 1 - 


K 


erf. 


±112 


exp(-a^a ) 


X cos(IIq') ((T^q; cos(ll(a) ± H sin(ll(a))) . (4-14) 

The equation for the order parameter 03.61) can be evaluated in the thermodynamic 
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limit to 


f = cos(f2a) exp 


0 ^, 0 ? 


(4.15) 


We have again omitted to write down the cumbersome expressions for the case of 
hnite TV, which nevertheless can readily be put into a numerical programme. 

Figure [13] shows a remarkable skill of the collective coordinate approach to describe 
the onset of global synchronisation and the order parameter f. The critical coupling 
strength for the global synchronisation at iVc = 1-70 is well captured. Furthermore, 
hnite-size scaling can be described within our framework as shown in Figure [T4| where 
we show a comparison of the critical coupling strength Kc{N) as calculated via our 
collective coordinate approach for variable network sizes N and results from direct 
simulations of the Kuramoto model fl2.ll) . As before we use as a proxy for the critical 
coupling strength the smallest value of the coupling strength K such that f > 0.8. 
The normalised size ^domain of the globally synchronised cluster, which we determine 
as the largest number of nodes for which non-trivial stationary solutions a exist, 
is depicted in Figure [151 The smooth gradual decrease of Tdomain with decreasing 
coupling strength iV, is replaced here by a different behaviour caused by the standing 
wave state and the partial synchronisation of oscillators with native frequencies close 
to ±11. 

Oscillators with native frequencies u ~ ±0 near the maxima of the native fre¬ 
quency distribution experience local synchronisation similar to the case of unimodally 
distributed native frequencies discussed in Section 14.21 In the case of a bimodal fre¬ 
quency distribution this leads to two partially synchronised clusters - one with fre¬ 
quency close to —O and another one with frequency close to ±0 (cf. Figure [T^ . 
With increasing coupling strength K the two clusters grow in size and will start to 
interact before, upon further increasing iV, they will merge at the onset of global 
synchronisation. We recall that this scenario only occurs provided the two peaks of 
the distribution of the native frequencies are sufficiently far separated allowing for a 
range in K for which they can partially synchronise without interacting too strongly 


12l | to form the standing wave state. We now set out to describe the standing wave 


state in our collective coordinate approach. 

In order to describe the effect of two partially synchronised clusters which rotate 
with non-uniform angular speeds of opposite direction we modify our ansatz and 
introduce a time-dependent phase function f{t) as an additional collective coordinate. 
We split the phase oscillators into two groups, one group (p~ describing the cluster 
centred around —hi, and one group (pf describing the cluster centred around ±12. We 
make the ansatz 


p>f = a{t){utT^)± fit) , (4.16) 

where uf are the native frequencies of the nodes participating in the cluster centred 
around ±12. Motivated by the results from direct simulations of the Kuramoto model 
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Figure 12: Snapshot of the phases ipj for K = 1.2 for the Kuramoto model with all- 
to-all coupling and bimodally distributed native frequencies with distribution fl4.13p . 
One can see clearly the two partially synchronised clusters with frequencies centred 
around O = ±0.75. The two clusters rotate with angular velocities of opposite sign 
forming a standing wave state. 

we assume that each of the clusters consists of N 2 < N/2 oscillators. Projecting the 
error onto the restricted subspace spanned by fl4.16p . i.e. onto dcpf/da = {ui =F O) 
and onto dipf/df = ±1, yields the desired evolution equations for a{t) and f{t). 
Projecting onto d(p~/da and dip~/df yields 

K 1 

« = 1 + ^) sin- ujr ) 

^ * j 

N2 N2 

- + ^) sin(a(u;7 ± wr) + 2afl - 2/) 

2 i j 

1 fC 

f = -^7 ) 

^ i,j 

1 K 

+ 2 S sin(a(a;" ± u;") ± 2afl - 2/) , 

^ i,j 


(4.17) 


(4.18) 


where here 


N 2 

+ ^)^ • (4-19) 

j 

The sums are taken over indices representing the nodes within the clusters ip~ (cf. 
fl4.10p h Due to symmetry projecting onto dipf /da and d'pj'/df reduces to the same 
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equation. The first sum in the right hand side of fl4.17p describes the interaction of 
oscillators within the partially synchronised cluster whereas the second sum describes 
the interaction of oscillators of one cluster with those of the respective other cluster. 

In the thermodynamic limit N ^ oo, the evolution equations for the collective 
coordinates simplify in the case when N 2 = N/2 to 

K 

a = 1 — — exp {—a^a^)a{l + cos(2/)) (4.20) 

f = sin(2/) . (4.21) 

Whereas in the case of global synchronisation the collective coordinate evolves to a 
stationary value, in the standing wave regime solutions of the system fl4.17p - fl4.18p 
or fl4.20p - fl4.2ip are oscillatory. These solutions can be found numerically. The order 
parameter can then be calculated as an average of flO.Op over one period Tp of the 
phase function f{t) and is given in the thermodynamic limit as 

f=- exp(-^) cos fit) dt , (4.22) 

Jo ^ 

In the thermodynamic limit the period Tp can be determined analytically. Dehning 
the collective coordinate a as an average over the period Tp, the Adler equation fl4.2ip 
can be solved analytically as 


f{t) = arctan 


/A + tan 

V ; 


(4.23) 


with A = iK/2) exp(—a^cr^). The associated period Tp is then dehned as 


T = 

-Lp 


df 


TT 


Vt- A sin(2/) yjVT - A^ 


(4.24) 


Note that because there are two counter-rotating clusters, the integration only goes 
to TT rather than to 27r. 

In Figure [T3] we show results of the collective coordinate approach for the order 
parameter f as a function of the coupling strength K. In practice we hrst test for 
global synchronisation, and if this cannot be achieved for any domain length ^domain, 
we test for the standing wave state. We have again allowed for local synchronisation 
whereby not all of the N/2 oscillators are synchronised (cf. Figure IT^ analogously 
to fl4.9p and fl4.10l) . The onset of the standing wave state at Kp = 1.05 is very well 
captured. The size of the synchronised clusters is shown in Figure [15] where we count 
the total sum of locally synchronised oscillators ip~ and in the case of the standing 
wave state for iF < 1.7. 

In Figure Uni we show a comparison of the actual temporal evolution of individual 
oscillators in the global synchronisation regime at iF = 2.5 and in Figure [TTl in the 
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Figure 13: Order parameter f as a function of the coupling strength K for a net¬ 
work with bimodal native frequency distribution. Depicted are results from a direct 
numerical integration of the Kuramoto model (12.Ih with N = 500 (continuous line, 
online red) and from the collective coordinate approach. We show results for global 
synchronisation (crosses, online blue) and for the standing wave state (open circles, 
online cyan). 


standing wave regime aX K = 1.1. The phases of the drifting oscillators which are 
not included within the collective coordinate analysis, are plotted simply by assuming 
that they are oscillating with their native frequencies. The actual phase dynamics of 
the synchronized oscillators is well described by our collective coordinate approach. 
One sees clearly the oscillatory behaviour of the phases in the standing wave regime 
which is caused by the interaction of the two counter-rotating clusters. The oscillation 
with period Tp = 5.8 is well captured by the dynamics of the collective coordinates 
and matches approximately the analytically obtained period Tp = 5.6 if we use the 
sample mean and variance of the native frequencies instead of D and in fl4.24p . 

5 Summary and Discussion 

The collective coordinate approach we propose allows for the description of networks 
of N oscillators. The dimension N is drastically reduced to a few n judiciously cho¬ 
sen collective coordinates; here we presented examples with n = 1 and n = 2. The 
approach is not restricted to the thermodynamic limit of inhnite network size and al¬ 
lows to study hnite networks. The approach can be used to study the synchronisation 
behaviour of networks, both global and partial, and determine the order parameter 
and the size of the synchronised clusters. Besides capturing this collective behaviour 
of oscillators the collective coordinate approach also is able to resolve the temporal 
evolution of individual oscillators for a wide range of coupling strength. 

We have corroborated our approach for the Kuramoto model with all-to-all cou¬ 
pling in numerical simulations for different distributions of the native frequencies. We 
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Figure 14: Scaling of the critical coupling strength for the onset of global syn¬ 
chronisation as a function of the network size N for a network with bimodal native 
frequency distribution. Depicted are results from a direct numerical integration of 
the Kuramoto model JED (crosses, online red) and from the collective coordinate 
approach (continuous line, online blue). The direct numerical simulations scale with 
a slope of 0.89. 
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Figure 15: Normalised length of phase synchronised domain as a function of the 
coupling strength for a network with bimodal native frequency distribution, calculated 
using the collective coordinate approach. The globally synchronised branch with 
-^domain = 1 Hs preceded for iF < 1.7 by a standing wave state in which, for K close 
to 1.7, all oscillators are involved (i.e. Tdomain = 1), but where the two partially 
synchronised clusters are not oscillating in phase. 
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Figure 16: Phases (p{t) calculated from simulations of the full Kuramoto model 
fl2.ll) (continuous lines or open circles, online red) and from the corresponding 1- 
dimensional system fl3.2p for the collective coordinate with ipi = a(t)ui (crosses, 
online blue) for a network of iV = 501 oscillators with a bimodal distribution fl4.13p 
of the native frequencies at coupling strength K = 2.5 corresponding to global syn¬ 
chronisation. Top: Temporal evolution of :p 47 o(t) for initial conditions v^j(O) = aoooi 
with ao = 0.5. The native frequency is a;47o = 1.11.Bottom: Snapshot of the phases 
(Pi{T) at time T = 20. 

found good agreement of our reduced 1-dimensional model (or 2-dimensional model 
in the case of bimodal native frequency distributions) with the full iV-dimensional 
system. In particular, the behaviour of the order parameter was well captured and 
the approach was able to describe soft second-order as well as explosive hrst-order 
transitions to synchronisation. We have illustrated that the collective coordinate ap¬ 
proach reproduces hnite size scalings of the full system. Furthermore, the approach 
allowed to describe the interplay between a standing wave state involving partially 
synchronised counter-rotating clusters and global synchronisation in networks with 
bimodal distribution of native frequencies. We have shown that the collective coor¬ 
dinates are able to capture the dynamics of individnal oscillators which is a much 
stronger form of approximation than just reproducing the collective behaviour. 

It is pertinent to caution that the method is by no means rigorous. The choice 
of collective coordinates is so far limited to a priori information obtained from direct 
numerical simulations of the full dynamical network. We have seen that transitory 
temporal evolution of oscillators in a Kuramoto model is only well described by the 
collective coordinate method provided the initial conditions are sufficiently close to 
the synchronisation manifold. Furthermore, the temporal evolution of individual os¬ 
cillators at the edge of a synchronised cluster is not accurately captured. To put our 
ansatz on a hrm theoretical footing which allows to describe its limitations is an open 
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Figure 17: Phases (p{t) calculated from simulations of the full Kuramoto model 
fl2.ll) (continuous lines or open circles, online red) and from the corresponding 2- 
dimensional system 04.211) for the collective coordinate with Lfi = ±/(t) 

(crosses, online blue) for a network of iV = 501 oscillators with a bimodal distribu¬ 
tion 04.131) of the native frequencies at coupling strength K = 1.3 corresponding to 
the standing wave regime. The size of the two respective counter-rotating clusters is 
N 2 = 210. Top: Temporal evolution of (^i8o(^) for initial conditions V5i(0) = ao(ci;i-|-r2) 
with ao = 0.5, i.e. a(0) = Oq and /(O) = 0. The native frequency is Wigo = —0.57. 
Bottom: Snapshot of the phases (pi{T) at time T = 20. 
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question. 


From a practical point of view, there are several issues which require further at¬ 
tention and which we plan to pursue in future research. First of all, whereas the 
general framework of collective coordinates is formulated for general network topolo¬ 
gies, we have only presented numerical results for the case of an all-to-all coupling. 
It is an interesting and important question to see whether the success of the method 
translates to more complex network topologies. 

Second, it is by no means clear that our ansatz captures all possible attractors of 
the full dynamical system. For exarnple, there are examples of networks where the 
Ott-Antonson method of reduction does not account for the actual dynamical 
behaviour observed in these networks (see the discussion in Martens et al. 0 ). In 
particular, chaotic dynamics is excluded from their framework. The collective coordi¬ 
nate approach is, in principle, capable of recovering chaotic dynamics by considering 
at least three collective coordinates. To test whether it actually is able to describe 
more complex dynamic behaviour is an interesting avenue to pursue. 

Thirdly, the success in describing the interaction between two partially synchronised 
clusters in the case of bimodally distributed native frequencies suggests that collec¬ 
tive coordinates may be used to reduce complex networks involving several clusters 
or communities. 

Fourthly, as we have seen in the numerical simulation, the collective coordinate ap¬ 
proach does not capture the interaction between the drifter oscillators and the syn¬ 
chronised oscillators. This leads to the collective coordinate behaviour not being able 
to accurately capture the oscillators which sit on the edge of the cluster. At a next 
step one can extend the approach to include drifters. 
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